** Below are programs built for all analysis used
set more off

** Calculate effect of tech on choice between vignettes **
cap prog drop domarginswithcurrentcomm
prog domarginswithcurrentcomm , rclass
	syntax [varlist(default=none)] [if] , [EXPression(string) dydx(string) at1(string) at2(string) statname(name) fmt(string) marginvar(string) post]
	if "`fmt'"!="" {
		assert "`statname'"!=""
		loc fmtstr fmt(`fmt')
	}
	if "`marginvar'"!="" {
		assert "`statname'"!="" | "`post'"=="post"
		loc marginvarstr marginvar(`marginvar')
	}
	if "`at1'"!="" {
		loc atstr at(`at1')
	}
	if "`at2'"!="" {
		assert "`at1'"!=""
		if "`marginvar'"!="" {
			disp "margin taken from at2; marginvar not used."
		}
		if "`dydx'"!="" {
			loc dydxeqname `dydx':
		}
		loc atstr `atstr' at(`at2')
	}
	preserve    // 88: add "restore" at end?
		qui replace comm = currentcomm if opt!=0
		disp "run this?"
		loc marginsstring margins `varlist' `if' , expression(`expression') dydx(`dydx') `atstr' `post'
		disp `"`marginsstring'"'
		`marginsstring'
		
 		matrix A0 = r(table)
		
		return add
end

** Logit marginal effects **
capture program drop mylogit_mle
program define mylogit_mle
args lnf Xb
	qui replace `lnf' = ln(invlogit(`Xb')*(1-$ML_y2) + (1-invlogit(`Xb'))*($ML_y2)) if $ML_y1==1
	qui replace `lnf' =  ln(invlogit(`Xb')*($ML_y2) + (1-invlogit(`Xb'))*(1-$ML_y2)) if $ML_y1==0
end	


** Report regression estimates **
capture prog drop svyregstat
program svyregstat
  local matrix A0 = r(table)
  matrix list A0
  local ls0 : display %5.2f A0[1,2] 
  local ls0p : display %5.3f A0[4,2]
  local lslc0 : display %5.2f A0[5,2]
  local lsuc0 : display %5.2f A0[6,2]
 
 nlcom cwd : -_b[1.tech] / _b[c.comm]
  matrix A = r(table)
  matrix list A
  local cwd : display %5.2f A[1,1] 
  local cwdp : display %5.3f A[4,1]
  local cwdlc : display %5.2f A[5,1]
  local cwduc : display %5.2f A[6,1]
  
  local subj : display %6.0f  e(N_clust)
  local ll : display %9.1f e(ll) 
  
 outreg , $outregopt ///
  ctitle("VARIABLES" , "Driving knowledge low") ///
  addrows("Subjects", "`subj'" \ "LogL", "`ll'" \ ///
   "CWD", "`cwd'" \ "", "[`cwdp']" \ ///
   "Work", "`ls0'" \ "", "[`ls0p']" ) ///
  merge(regskill)  
/* Syntax  
  1 = file handle for graph data
  2 = skill
  3 = column title
  4 = file handle for outreg */
end


** Report regression estimates by skill **
capture prog drop svyregstatskill
program svyregstatskill
  local matrix A0 = r(table)
  matrix list A0
  local ls0 : display %5.2f A0[1,3] 
  local ls0p : display %5.3f A0[4,3]
  local ls0lc : display %5.2f A0[5,3]
  local ls0uc : display %5.2f A0[6,3]
  local ls1 : display %5.2f A0[1,4] 
  local ls1p : display %5.3f A0[4,4]
  local ls1lc : display %5.2f A0[5,4]
  local ls1uc : display %5.2f A0[6,4]
  
	
 nlcom cwd0 : -_b[1.tech] / _b[c.comm]
  matrix A0 = r(table)
  matrix list A0
  local cwd0 : display %5.2f A0[1,1] 
  local cwd0p : display %5.3f A0[4,1]
  local cwd0lc : display %5.2f A0[5,1]
  local cwd0uc : display %5.2f A0[6,1]

	
 nlcom cwd1 : - ( _b[1.tech] + _b[1.`2'#1.tech] ) / ///
   ( _b[c.comm] + _b[1.`2'#c.comm] )
  matrix A1 = r(table)
  matrix list A1
  local cwd1 : display %5.2f A1[1,1] 
  local cwd1p : display %5.3f A1[4,1]  
  local cwd1lc : display %5.2f A1[5,1]
  local cwd1uc : display %5.2f A1[6,1]


  local subj : display %6.0f  e(N_clust)
  local ll: display %9.1f e(ll) 
  
outreg , $outregopt ///
  ctitle("VARIABLES" , `3') ///
  addrows( "Subjects", "`subj'" \ "LogL", "`ll'" \ ///
   "CWD low skill", "`cwd0'" \ "", "[`cwd0p']" \ ///
   "CWD high skill", "`cwd1'" \ "", "[`cwd1p']" \ ///
   "Work low skill", "`ls0'" \ "", "[`ls0p']" \ ///
   "Work high skill", "`ls1'" \ "", "[`ls1p']" ) ///
  merge(`4')  
  
/* Syntax  
  1 = file handle for graph data
  2 = skill
  3 = column title
  4 = file handle for outreg */
end

   
** add regression estimates ** 
cap prog drop addvarstat
prog addvarstat
	syntax varname [if] , statname(name) [scale(real 1.)]
	marksample touse
	sum `varlist' if `touse' & e(sample) , meanonly
	estadd scal `statname' = r(mean) * `scale'
end

cap prog drop domarginswithcurrentcomm_old
prog domarginswithcurrentcomm_old , rclass
	syntax [varlist(default=none)] [if] , [EXPression(string) dydx(string) at(string) ESTADDstring(string asis)]
	preserve
		qui replace comm = currentcomm if opt!=0
		margins `varlist' `if' , expression(`expression') dydx(`dydx') at(`at')
		`estaddstring'
		return add
	restore
end

** rename regression estimates ** 
cap prog drop renameestimates
prog renameestimates , eclass
	syntax anything , newname(name) [modlocal(name) newlocalname(name) all]
	mat eb = e(b)
	loc cn: colnames eb
	loc cn: subinstr loc cn "`anything'" "`newname'" , `all'
	mat colnames eb = `cn'
	ereturn repost b=eb , rename
	if "`modlocal'"!="" {
		ereturn local `newlocalname' `e(`modlocal')'
	}
end

** add regression estimates ** 
cap prog drop adddeltamethodstat
prog adddeltamethodstat
	syntax , type(string) statname(name) [fmt(string) marginvar(string)]
	assert inlist("`type'","nlcom","margin","margins","lincom")
	if "`fmt'"=="" {
		loc fmt %9.2f
	}
	if "`marginvar'"=="" & inlist("`type'","margin","margins") {
		loc marginvar _cons
	}
	if "`type'"=="nlcom" {
		estadd scal `statname'_b = r(b)[1,1]
		estadd scal `statname'_se = sqrt( r(V)[1,1] )
		estadd scal `statname'_p = 2 * ttail( 1e9 , abs(e(`statname'_b)/e(`statname'_se)) )
	}
	else if "`type'"=="lincom" {
		estadd scal `statname'_b = r(estimate)
		estadd scal `statname'_se = r(se)
		estadd scal `statname'_p = r(p)
	}
	else if inlist("`type'","margin","margins") {
		estadd scal `statname'_b = r(table)["b","`marginvar'"]
		estadd scal `statname'_se = r(table)["se","`marginvar'"]
		estadd scal `statname'_p = r(table)["pvalue","`marginvar'"]	  
	}
	loc stars ""
	foreach threshold of numlist .1 .05 .01 {
		if e(`statname'_p)<=`threshold' {
			loc stars `stars'*
		}
	}
	estadd loc `statname' `=string( e(`statname'_b) , "`fmt'" )'`stars'		
end

cap prog drop donlcomandadd
prog donlcomandadd
	syntax , expression(string) [statname(name) fmt(string) marginvar(string)]
	if "`fmt'"!="" {
		assert "`statname'"!=""
		loc fmtstr fmt(`fmt')
	}
	if "`marginvar'"!="" {
		assert "`statname'"!=""
		loc marginvarstr marginvar(`marginvar')
	}
	nlcom `expression'
	if "`statname'"!="" {
		adddeltamethodstat , type(nlcom) statname(`statname') `fmtstr' `marginvarstr'
	}
end


** Predict inattention ** 
cap prog drop getinattentionprediction
prog getinattentionprediction
	syntax varlist(min=2) [if] , GENerate(name)
	marksample touse
	reg `varlist' if `touse' & inattcheckqn & opt==0
	predict `generate'
	sum `generate'
	if `r(min)'<0 {
		disp as err "Predicted inattention below zero, set to zero."
		replace `generate' = 0 if `generate'<0
	}
	if `r(max)'>1 {
		disp as err "Predicted inattention above one, set to one."
		replace `generate' = 1 if `generate'>1
	}
end

cap prog drop formatp
program define formatp
    syntax name
    
    local pvar `"`1'"'
    
    if (e(`pvar') < 0.001) {
        estadd loc `pvar' "<0.001", replace
    }
    else {
        local formatted_p = string(e(`pvar'), "%9.3f")
        estadd loc `pvar' "`formatted_p'", replace
    }
end

cap prog drop formatll
program define formatll
    scalar ll = e(ll)
    local ll_str = string(abs(ll), "%15.10f")
    local dot_pos = strpos("`ll_str'", ".")
	local first_nonzero = 0
	local len_after = strlen("`ll_str'") - `dot_pos'
	local after_dot = substr("`ll_str'", `dot_pos' + 1, .)
    forvalues i = 1(1)`len_after' {
        local char = substr("`after_dot'", `i', 1)
        if `char' != 0 {
            local first_nonzero = `i'
            continue, break
        }
    }
    if substr("`ll_str'", `dot_pos' + 1, 1) == "0" {
		estadd loc loglikelihood = string(ll, "%9.`first_nonzero'f"), replace
	}
	else {
		estadd loc loglikelihood = string(ll, "%9.2f"), replace
	}
end	


** logit estimates ** 
cap prog drop mylogit
prog mylogit
	syntax varlist(min=2 fv) [if] , ERRorvar(varname)
	gettoken y rhs : varlist
	marksample touse
	ml model lf mylogit_mle (`y' `errorvar' = `rhs') if `touse' , vce(robust)
	ml max, iterate(100)
end


capture program drop mylogit_mle
program define mylogit_mle
args lnf Xb
	qui replace `lnf' = ln(invlogit(`Xb')*(1-$ML_y2) + (1-invlogit(`Xb'))*($ML_y2)) if $ML_y1==1
	qui replace `lnf' =  ln(invlogit(`Xb')*($ML_y2) + (1-invlogit(`Xb'))*(1-$ML_y2)) if $ML_y1==0
end	


** analysing hrv data ** 

capture program drop hrv
program define hrv
	global ecg rmssd hrvti hr sdnn si lf_fft hf_fft ratio_fft lf_ar hf_ar ratio_ar sampen
	foreach hrv in $ecg {
		rename `hrv'`1' `hrv'3
		rename `hrv'`2' `hrv'4
		drop `hrv'*_*_*
		gen `hrv'_delta2 = (`hrv'4 - `hrv'3) / `hrv'3
	}
	label var rmssd_delta2 "Successive differences (root mean square)"
	label var hr_delta2 "Heart rate"
	label var sdnn_delta2 "Normal-to-normal (std deviation)"
	label var hrvti_delta2 "Heart Rate Variability Triangular Index"
	label var si_delta2 "Stress index (Baevsky's Stress Index)"
	label var lf_fft_delta2 "Low frequency power (Fast Fourier Transform)"
	label var hf_fft_delta2 "High frequency power (Fast Fourier Transform)"
	label var ratio_fft_delta2 "LF/HF ratio (Fast Fourier Transform)"
	label var lf_ar_delta2 "Low frequency power (autoregressive)"
	label var hf_ar_delta2 "High frequency power (autoregressive)"
	label var ratio_ar_delta2 "LF/HF ratio (autoregressive)"
	label var sampen_delta2 "Sample entropy"
	
	if word(subinstr("`1'", "_", " ", .), 1) == "1" {
		local ending1 DrvingInstructn
	}
	else if word(subinstr("`1'", "_", " ", .), 1) == "2" {
		local ending1 lsQns
	}
	else if word(subinstr("`1'", "_", " ", .), 1) == "3" {
		local ending1 RevlDestn
	}
	
	if word(subinstr("`1'", "_", " ", .), 2) != "0" {
		local sec12 "`=word(subinstr("`1'", "_", " ", .), 2)'s post-"
	} 
	if word(subinstr("`1'", "_", " ", .), 3) != "0" {
		local sec13 "`=word(subinstr("`1'", "_", " ", .), 3)'s pre-"
	} 
	
	if word(subinstr("`2'", "_", " ", .), 2) != "0" {
		local sec22 "`=word(subinstr("`2'", "_", " ", .), 2)'s post-"
	} 
	if word(subinstr("`2'", "_", " ", .), 3) != "0" {
		local sec23 "`=word(subinstr("`2'", "_", " ", .), 3)'s pre-"
	} 
	
	label var hr4 `"Heart rate (driving) `sec22'driving to `sec23'end"'
	label var rmssd4 "RMSSD (driving) `sec22'driving to `sec23'end"
	label var hf_fft4 "HF (driving) `sec22'driving to `sec23'end"
	label var ratio_fft4 "LF/HF ratio (driving) `sec22'driving to `sec23'end"

	label var hr3 `"Heart rate (pre-drive) `sec12'start to `sec13'`ending1'"'
	label var rmssd3 "RMSSD (pre-drive) `sec12'start to `sec13'`ending1'"
	label var hf_fft3 "HF (pre-drive) `sec12'start to `sec13'`ending1'"
	label var ratio_fft3 "LF/HF ratio (pre-drive) `sec12'start to `sec13'`ending1'"
	
	foreach x in hr4 rmssd4 hf_fft4 ratio_fft4 hr3 rmssd3 hf_fft3 ratio_fft3 stai {
	  local varlab : var label `x'
	  gen `x'_ln = ln(`x')
	  label var `x'_ln "`varlab' (ln)" 
	}
end

** Programs for field experiemnt ** 

capture program drop fieldregstat
program define fieldregstat
  local obs = e(N)
  local lnl : di %7.2f e(ll)
  local r2 : di %04.3f e(r2)
  unique id if e(sample)
  local subject : display %6.0f `r(unique)'
  summ `1' if e(sample) & tech == 0
  local dvmean : display %04.3f `r(mean)'
  
  margins, dydx(tech)  
  matrix A = r(table)
  matrix list A
  local margeff : display %5.3f A[1,1] 
  local pvalue : display %5.3f A[4,1]
 
  outreg , keep( `2' ) ctitle("VARIABLES" , `3' ) landscape ///
   plain coljust(lc) varlabel stats(b se p) sdec(3) ///
   nostar noautosumm ///
   addrows( "Observations" , "`obs'" \ "lnL" , "`lnl'" \ ///
    "R2" , "`r2'" \ "Subjects" , "`subject'" \ "All subjects" , " " \ ///
	"- Mean `1' (no tech)" , "`dvmean'" \ ///
	"- Marginal effect of tech" , "`margeff'" \ " " , "[`pvalue']" ) ///
	merge(`4')
end 

/* Syntax
   1 = dependent variable
   2 = explanatory variables to report
   3 = column title
   4 = name of outreg handle
*/
 
 
capture program drop fieldregstatskill
program define fieldregstatskill

  local obs = e(N)
  local lnl : di %7.2f e(ll)
  local r2 : di %04.3f e(r2)
  unique id if e(sample)
  local subject : display %6.0f `r(unique)'
  summ `1' if e(sample) & tech == 0 
  local dvmean : display %04.3f `r(mean)'
  summ `1' if e(sample) & selfseghigh == 0 & tech == 0 
  local dvmean0 : display %04.3f `r(mean)'
  summ `1' if e(sample) & selfseghigh == 1 & tech == 0 
  local dvmean1 : display %04.3f `r(mean)'
  
  margins, dydx(tech) at(selfseghigh == (0 1))
  matrix A = r(table)
  matrix list A
  local margeff0 : display %5.3f A[1,3] 
  local pvalue0 : display %5.3f A[4,3]
  local ls0lc : display %5.3f A[5,3] 
  local ls0uc : display %5.3f A[6,3]
  local margeff1 : display %5.3f A[1,4] 
  local pvalue1 : display %5.3f A[4,4]
  local ls1lc : display %5.3f A[5,4] 
  local ls1uc : display %5.3f A[6,4]
  
  file open regworkci using `2' , ///
    write replace  // compile data for Figure .. (next bookmark)
  file write regworkci "selfseghigh" _tab "ls" ///
    _tab "lc" _tab "uc" _newline
  file write regworkci "0" _tab "`margeff0'" _tab "`ls0lc'" ///
    _tab "`ls0uc'" _newline 
  file write regworkci "1" _tab "`margeff1'" _tab "`ls1lc'" ///
    _tab "`ls1uc'" _newline 
  file close regworkci

  if "`6'" == "" {
  local meanlbl `1'
  }
  else {
  local meanlbl "`6'"
  }

  outreg , keep( `3' ) ctitle("VARIABLES" , `4' ) landscape ///
   plain coljust(lc) varlabel stats(b se p) sdec(3) ///
   nostar noautosumm ///
   addrows( "Observations" , "`obs'" \ "lnL" , "`lnl'" \ ///
    "R2" , "`r2'" \ "Subjects" , "`subject'" \ ///
	"Low-skill" , " " \ ///
    "- low: Mean work (no tech)" , "`dvmean0'" \ ///
	"- low: Marg. effect of tech" , "`margeff0'" \ ///
    "  low: p-value" , "[`pvalue0']" \ ///
	"High-skill" , " " \ ///
    "- high: Mean work (no tech)" , "`dvmean1'" \ ///
	"- high: Marg. effect of tech" , "`margeff1'" \ ///
    "  high: p-value" , "[`pvalue1']" ) merge(`5')  
	
end 

/* Syntax
   1 = dependent variable
   2 = name of file to store coefficients
   3 = explanatory variables to report
   4 = column title
   5 = name of outreg handle
   6 = outcome (if dependent variable, omit) */

   capture program drop fieldregstatfe
program define fieldregstatfe
  local obs = e(N)
  local lnl : di %7.2f e(ll)
  local r2 : di %04.3f e(r2)
  unique id if e(sample)
  local subject : display %6.0f `r(unique)'
  summ `1' if e(sample) & tech == 0
  local dvmean : display %04.3f `r(mean)'
  
  margins, dydx(tech)  
  matrix A = r(table)
  matrix list A
  local margeff : display %5.3f A[1,1] 
  local pvalue : display %5.3f A[4,1]
 
  outreg , keep( `2' ) ctitle("VARIABLES" , `3' ) landscape ///
   plain coljust(lc) varlabel stats(b se p) sdec(3) ///
   nostar noautosumm ///
   addrows( "Group fixed effects" , `4' \ ///
    "Destination f.e." , `5' \ ///
    "Destination x skill f.e." , `6' \ ///
    "Observations" , "`obs'" \ "lnL" , "`lnl'" \ "R2" , "`r2'" \ ///
    "Subjects" , "`subject'" \ "All subjects" , " " \ ///
	"- Mean `1' (no tech)" , "`dvmean'" \ ///
	"- Marginal effect of tech" , "`margeff'" \ " " , "[`pvalue']" ) ///
	merge(`7')
end 

/* Syntax
   1 = dependent variable
   2 = explanatory variables to report
   3 = column title
   4 = whether included group fixed effects
   5 = whether included destination fixed effects
   6 = whether included destination x skill fixed effects
   7 = name of outreg handle */
   
capture program drop fieldregstatskillfe
program define fieldregstatskillfe

  local obs = e(N)
  local lnl : di %7.2f e(ll)
  local r2 : di %04.3f e(r2)
  unique id if e(sample)
  local subject : display %6.0f `r(unique)'
  summ `1' if e(sample) & tech == 0 
  local dvmean : display %04.3f `r(mean)'
  summ `1' if e(sample) & selfseghigh == 0 & tech == 0 
  local dvmean0 : display %04.3f `r(mean)'
  summ `1' if e(sample) & selfseghigh == 1 & tech == 0 
  local dvmean1 : display %04.3f `r(mean)'
  local modeltype e(cmd)
  

	if inlist(`modeltype', "regress", "reghdfe"){
		lincom 1.tech
		local margeff0 : display %5.3f r(estimate)
		local pvalue0  : display %5.3f r(p)
		local ls0lc    : display %5.3f r(lb)
		local ls0uc    : display %5.3f r(ub)

		lincom 1.tech + 1.tech#1.selfseghigh
		local margeff1 : display %5.3f r(estimate)
		local pvalue1  : display %5.3f r(p)
		local ls1lc    : display %5.3f r(lb)
		local ls1uc    : display %5.3f r(ub)
	} 
	else {
		margins, dydx(tech) at(selfseghigh = (0 1))
		matrix A = r(table)
		matrix list A
		local margeff0 : display %5.3f A[1,3]
		local pvalue0  : display %5.3f A[4,3]
		local ls0lc    : display %5.3f A[5,3]
		local ls0uc    : display %5.3f A[6,3]
		local margeff1 : display %5.3f A[1,4]
		local pvalue1  : display %5.3f A[4,4]
		local ls1lc    : display %5.3f A[5,4]
		local ls1uc    : display %5.3f A[6,4]
	}

  
  file open regworkci using `2' , ///
    write replace  // compile data for Figure .. (next bookmark)
  file write regworkci "selfseghigh" _tab "ls" ///
    _tab "lc" _tab "uc" _newline
  file write regworkci "0" _tab "`margeff0'" _tab "`ls0lc'" ///
    _tab "`ls0uc'" _newline 
  file write regworkci "1" _tab "`margeff1'" _tab "`ls1lc'" ///
    _tab "`ls1uc'" _newline 
  file close regworkci

  if "`9'" == "" {
  local meanlbl `1'
  }
  else {
  local meanlbl "`9'"
  }
  
  outreg , keep( `3' ) ctitle("VARIABLES" , `4' ) landscape ///
   plain coljust(lc) varlabel stats(b se p) sdec(3) ///
   nostar noautosumm ///
   addrows( "Group fixed effects" , `5' \ ///
    "Destination f.e." , `6' \ ///
    "Destination x skill f.e." , `7' \ ///
    "Observations" , "`obs'" \ "lnL" , "`lnl'" \ "R2" , "`r2'" \ ///
    "Subjects" , "`subject'" \ ///
	"Low-skill" , " " \ ///
    "- low: Mean work (no tech)" , "`dvmean0'" \ ///
	"- low: Marg. effect of tech" , "`margeff0'" \ ///
    "  low: p-value" , "[`pvalue0']" \ ///
	"High-skill" , " " \ ///
    "- high: Mean work (no tech)" , "`dvmean1'" \ ///
	"- high: Marg. effect of tech" , "`margeff1'" \ ///
    "  high: p-value" , "[`pvalue1']" ) merge(`8')  
end 

/* Syntax
   1 = dependent variable
   2 = name of file to store coefficients
   3 = explanatory variables to report
   4 = column title
   5 = whether included group fixed effects
   6 = whether included destination fixed effects
   7 = whether included destination x skill fixed effects
   8 = name of outreg handle
   9 = outcome (if dependent variable, omit) 
*/



capture program drop fieldmedianstat
program define fieldmedianstat
  local obs = e(N)
  local lnl : di %7.2f e(ll)
  local r2 : di %04.3f e(r2_p)
  unique id if e(sample)
  local subject : display %6.0f `r(unique)'
  summ `1' if e(sample) & tech == 0 , d
  local dvmedian : display %04.3f `r(p50)'
  
  margins, dydx(tech)  
  matrix A = r(table)
  matrix list A
  local margeff : display %5.3f A[1,1] 
  local pvalue : display %5.3f A[4,1]
 
  outreg , keep( `2' ) ctitle("VARIABLES" , `3' ) landscape ///
   plain coljust(lc) varlabel stats(b se p) sdec(3) ///
   nostar noautosumm ///
   addrows( "Observations" , "`obs'" \ "lnL" , "`lnl'" \ ///
    "Pseudo R2" , "`r2'" \ "Subjects" , "`subject'" \ "All subjects" , " " \ ///
	"- Median `1' (no tech)" , "`dvmedian'" \ ///
	"- Marg. effect of tech" , "`margeff'" \ " " , "[`pvalue']" ) ///
	merge(`4')
end 

/* Syntax
   1 = dependent variable
   2 = explanatory variables to report
   3 = column title
   4 = name of outreg handle
*/
 
